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ABSTRACT 

In this work we present for the first time an analytic framework for calculating the individ- 
ual and joint distributions of the n-th most massive or n-th highest redshift galaxy cluster 
for a given survey characteristic allowing to formulate ACDM exclusion criteria. We show 
r"! i that the cumulative distribution functions steepen with increasing order, giving them a higher 

constraining power with respect to the extreme value statistics. Additionally, we find that the 
order statistics in mass (being dominated by clusters at lower redshifts) is sensitive to the 
matter density and the normalisation of the matter fluctuations, whereas the order statistics in 
^ redshift is particularly sensitive to the geometric evolution of the Universe. For a fixed cos- 

^ I mology, both order statistics are efficient probes of the functional shape of the mass function 

at the high mass end. To allow a quick assessment of both order statistics, we provide fits as a 
function of the survey area that allow percentile estimation with an accuracy better than two 
per cent. Furthermore, we discuss the joint distributions in the two-dimensional case and find 
that for the combination of the largest and the second largest observation, it is most likely to 
find them to be realised with similar values with a broadly peaked distribution. When combin- 
^ , ing the largest observation with higher orders, it is more likely to find a larger gap between the 

\^ ■ observations and when combining higher orders in general, the joint pdf peaks more strongly. 

Having introduced the theory, we apply the order statistical analysis to the SPT massive 
cluster sample and MCXC catalogue and find that the ten most massive clusters in the sample 
are consistent with ACDM and the Tinker mass function. For the order statistics in redshift, we 
find a discrepancy between the data and the theoretical distributions, which could in principle 
indicate a deviation from the standard cosmology. However, we attribute this deviation to the 
uncertainty in the modelling of the SPT survey selection function. In turn, by assuming the 
ACDM reference cosmology, order statistics can also be utilised for consistency checks of the 
completeness of the observed sample and of the modelling of the survey selection function. 
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1 INTRODUCTION 

Clusters of galaxies represent the top of the hierarchy of gravita- 
tionally bound structures in the Universe and can be considered 
as tracers of the rarest peaks of the initial density field. This fea- 
ture renders their abundance across the cosmic history a valuable 
probe of cosmology (for an overview of cluster cosmology see 
e.g. Voit 2005; Allen et al. 2011, and references therein). The re- 
cent years brought significant advances to the field from an ob- 
servational point of view. Past and present surveys, like e.g. the 
ROSAT All Sky Survey (RASS; Voges etal. 1999), the Massive 
Cluster Survey (MACS; Ebeling et al. 2001) and the Southpole 
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Telescope (SPT; Carlstrom et al. 2011), provided rich data for a 
multitude of massive clusters (> 10^^ M©). In the near future, 
cluster data will be drastically extended in terms of complete- 
ness, coverage and depth by surveys like for instance PLANCK 
(Tauber, J. A. et al. 2010), eROSITA (Cappelluti et al. 2011) and 
EUCLID (Laureijs et al. 2011), allowing for statistical analyses of 
the samples with increasing quality. 

A particular form of statistical analysis that recently en- 
tered focus are falsification experiments of the concordance 
ACDM cosmology, based on the discovery of a single (or a 
number) of cluster(s) being so massive that it (they) could 
not have formed in the standard picture (Hotchkiss 2011; 
Hoyle et al. 201 1 ; Mortonson et al. 201 1 ; Harrison & Coles 
2012; Harrison & Hotchkiss 2012; Holz & Perlmutter 2012; 
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Waizmann et al. 2012a,b). These studies were triggered by the 
discovery of massive clusters at high redshift (see e.g. MuUis et al. 
2005; Jeeetal. 2009; Rosati et al. 2009; Foley et al. 2011; 
Menanteau et al. 2012; Stalder et al. 2012). 

However, the usage of a single observation for such falsifi- 
cation experiments requires statistical care since several subtleties 
have to be taken into account. From the theoretical point of view, 
it is necessary to include the Eddington bias (Eddington 1913) in 
mass, as discussed in Mortonson et al. (2011) and the bias that 
stems from the a posteriori choice of the redshift interval for the 
analysis (Hotchkiss 2011). From the observational point of view, it 
might, particularly for very high redshift systems, be difficult to de- 
fine the survey area and selection function that are appropriate for 
the statistical analysis. Combining all of these effects, recent studies 
(Hotchkiss 2011; Harrison & Coles 2012; Harrison & Hotchkiss 
2012; Waizmann et al. 2012a,b) converge to the finding that, when 
taken alone, none of the single most massive known clusters can be 
considered in tension with the concordance ACDM cosmology. 

Conceptually, inference based on a single observation is not 
desirable, because by nature the extreme value might not be rep- 
resentative for the underlying distribution from which it is suppos- 
edly drawn. Thus, it is advised to incorporate statistical information 
from the sample of the most massive high redshift clusters, which 
in turn are also particularly sensitive to the underlying cosmologi- 
cal model since they probe the exponentially suppressed tail of the 
mass function. 

In this work, we introduce order statistics as a tool for analyt- 
ically deriving distribution functions for all members of the mass 
and redshift hierarchy ordered by magnitude. By dividing our anal- 
ysis in the observables mass and redshift, we avoid the bias due to 
an a posteriori definition of redshift intervals (Hotchkiss 201 1) and 
avoid as well the arbitrariness of an a priori choice that had been 
necessary in our previous works based on the extreme value statis- 
tics. Furthermore, the formalism also allows for the formulation of 
joint probabilities of the order statistics. In the second part of this 
work, we compare our individual and joint analytic distributions to 
observed samples of massive galaxy clusters. 

This paper is structured according to the following scheme. 
In Sect. 2, we introduce the statistical branch of order statistics by 
discussing the basic mathematical relations in Sect. 2.1 and by ap- 
plying the formalism to the distribution of massive galaxy clusters 
in mass and redshift in Sect. 2.2. This is followed by a discussion 
of how the order statistics of haloes in mass and redshift depends 
on cosmological parameters in Sect. 3. In order to compare our an- 
alytic results to observations, we prepare observed cluster samples 
for the analysis in Sect. 4. Afterwards, we discuss the results of the 
comparison for the case of the individual order statistic in Sect. 5 
and for the joint case in Sect. 6. Then, we summarise our findings 
in Sect. 7 and draw our conclusions in Sect. 8. In Appendix A we 
give a more detailed overview of order statistics and in Appendix B 
fitting formulae for the order statistics in mass and redshift are pre- 
sented. 

Throughout this work, unless stated otherwise, we adopt the 
Wilkinson Microwave Anisotropy Probe 7-year (WMAP7) param- 
eters {a^o,a^o,a^o,h,o-^) = (0.727,0.273,0.0455,0.704,0.811) 

(Komatsu etal. 2011). 



2 ORDER STATISTICS 

Order statistics (for an introduction, see e.g. Arnold et al. 1992; 
David & Nagaraja 2003) is the study of the statistics of ordered 



(sorted by magnitude) random variates. In this section, the basic 
mathematical relations and the connection to cosmology are intro- 
duced as they will be needed in remainder of this work. 

2.1 Mathematical prerequisites 

Let Xi, X2, . . . , be a random sample of a continuous population 
with the probability density function (pdf), f(x), and the corre- 
sponding cumulative distribution function (cdf), F(x). Further, let 
< X(2) < • • • < be the order statistic, the random variates 
ordered by magnitude, where is the smallest (minimum) and 
X(n) denotes the largest (maximum) variate. It can be shown (see 
Sect. Al) that the pdf of (1 < / < ^7) is given by 



■ [F(x)r'[i-F(x)r-^f(x). 



(i-\)l(n-i)l 
The corresponding cdf of the i-ih order reads then 



(1) 



(2) 



and the distribution function of the smallest and the largest value 
are found to be 



F(i)(x) = 1 - [1 - Fix)]" , 
and 

FUx) = [Fix)]" . 



(3) 



(4) 



In the limit of very large sample sizes both F(„)(x) and F(i)(x) can 
be described by a member of the general extreme value (GEV) dis- 
tribution (Fisher & Tippett 1928; Gnedenko 1943) 



G(x) - exp 



I + y 



-i/r 



(5) 



where a is the location-, /3 the scale- and 7 is the shape-parameter. 
Usually these parameters are obtained directly from the data or 
from an underlying model (see for instance Coles (2001)). 

Apart from the distributions of the single order statistics, it 
is very interesting to derive joint distribution functions for several 
orders. The joint pdf of the two order statistics X(r),X(s) (1 < r < 
s ^ n) is for x < y given by (see Appendix A for a more detailed 
discussion) 

n\ 

f{r){s)ix,y) -- 



{r - l)\{s - r - l)\{n - s)\ 
X lF(x)Y-' [F(y) - F(x)y-'-' [1 - Fiy)^ 
xmfiy). 



(6) 



The joint cumulative distribution function can e.g. be obtained by 
integrating the pdf above or by a direct argument and is found to be 
given by 



(;• -(•)!(«-;■)! 

j=s i=r ^ ^ ' 

X \F{x)\ {Fiy) - F{x)f [1 " ^(j)]"'' • (7) 

Analogously the above relations can be generalised to the joint pdf 
of , . . . , Xn], (1 ^ '^i < • • • < n]^ ^n)iox x\ ^ • • • ^ Xk, which is 
given by 

n\ 



(ni - l)\(n2 -ni -l)l---(n- nk)\ 

X [F{x,)r-' fix,) [F{X2) - F{x,)r-''^-' 

x/(x2)---[l-F(x,)f-"V(x,). (8) 
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Figure 1. Cumulative distribution functions of the first fifty orders from F(^n-A9) to F{n) in mass (left panel) and in redshift (right panel). For comparison, the 
GEV distribution of the maxima is shown by the red, dashed line for both cases. All distributions were calculated for the full sky, assuming the Tinker mass 
function. For the order statistics in redshift, a limiting mass of mum = 10^^ M© has been adopted. 



Further details and derivations concerning order statistics can be 
found in the Appendix A. In the remainder of this work we will 
repeatedly make use of percentiles. In statistics, a percentile is de- 
fined as the value of a variable below which a certain percentage, /?, 
of observations fall. Percentiles can be directly obtained from the 
inverse of the cdf and will be hereafter denoted as Qp. 



2.2 Connection to cosmology 

As outlined in the previous subsection, the only quantity that is 
needed for calculating the cdfs, F(j^{x), of the order statistics (see 
equation 2) is the cdf, F{x), of the underlying distribution from 
which the sample is drawn. Assuming the random variates, X/, to be 
the masses of galaxy clusters, then the cdf, F(m), can be calculated 
(see e.g. Harrison & Coles (2012)) by means of 



F{m) = 



/s 



sky 



Jo 



dzdM 



dV n{M,z) 



dz dM 



where the total number of clusters, Mot, is given by 
dV n{M,z) 



Mot = /s 



sky 



Jo Jo 



dzdM 



dz dM 



(9) 



(10) 



Here, /sky is the fraction of the full sky that is observed, (dV/dz) 
is the volume element and n(m,z) is the halo mass function. If 
needed, the corresponding pdf can always be obtained by f(m) = 
dF(m)/dm. 

Analogously, the order statistics can be calculated as well for 
the redshift instead of the mass. In this case the cdf reads 

^sky 



F(z)- 



Jo Jm 



dzdM 



dV n{M,z) 



dM 



(11) 



where 



■fs 



sky 



Iff 



dzdM 



dV n{M,z) 



dz dM 



(12) 



For the latter, the order statistics does no longer depend only on 
the survey area via /sky but, in addition the selection function of the 
survey has to be included via a limiting survey mass, miim(z). In this 
work, we do not attempt to model the possible redshift dependence 
of miim and assume it to be constant throughout the remainder of 
this work. 



With the distributions F{m) and F{z) at hand, we can now eas- 
ily derive the cdfs of the corresponding order statistics. Since we 
will focus in this work on the few largest values, we will refer to 
the distribution of the maximum, F{n){x), as first order, to the sec- 
ond largest as second order and so on. 

We calculated the distributions of the first fifty orders from 
F(„)(x) to F(„_49)(x), where n - Mot, and present the results in Fig. 1 
for the mass (left panel) and redshift (right panel). In both panels 
the color decodes the order of the distribution, ranging from the 
blue for F(n){x) to the green for F(n_^^){x). For both cases we as- 
sumed /sky = 1, a redshift range of < z < oo and the Tinker et al. 
(2008) mass function. In the case of the order statistics in redshift, 
we assume a limiting survey mass of = 10^^ M©. It can be 
nicely seen how, with increasing order (from blue to green), the 
cdfs shift in both cases to smaller values of the mass or redshift. 

A first important result is that, with the increasing order, the 
cdfs steepen, which results in an enhanced constraining power, 
since small shifts in the mass or redshift may yield large differences 
in the derived probabilities. In this sense the higher orders will be 
more useful for falsification experiments than the extreme value 
distribution which, due to its shallow shape, requires extremely 
large values of the observable to statistically rule out the underly- 
ing assumptions. Since higher orders encode information from the 
n most extreme objects, deviations from the expectation are sta- 
tistically more significant for n values instead of a single extreme 
one. 

In addition, we compare the distribution of the maxima 
F(„)(m) and F(„)(z) to those obtained from a extreme value ap- 
proach (Davis et al. 201 1 ; Waizmann et al. 2012a, Metcalf & Waiz- 
mann in prep.) based on the void probability (White 1979), using 
equation 5. For both cases presented in Fig. 1, the red, dashed curve 
of the GEV distribution, G(x), agrees very well with the directly 
calculated F(^n){^)- 

In order to allow a quick estimation of the distributions of the 
order statistics, we provide in the Appendix B also fitting formulae 
for F{x) as a function of the survey area for the cases of mass and 
redshift. The fitting formulae for the distribution in mass allow an 
estimation of the quantiles in the range from the 2-percentile, 22, to 
the 98-percentiles, g98, with an accuracy better than one per cent 
for As > 200 deg^ and for the ten largest masses. In the instance of 
the order statistics in redshift, the quality of the fits depends on 
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as well. For = 10^^ M© an accuracy of better than two per cent 
can be achieved for As > 2000 deg^ and for mum = 5 x 10^"^ M© 
the same accuracy is obtained down to As = 100 deg^. A more 
detailed discussion of the fitting functions and their performance 
can be found in Appendix B. 

In the remaining part of this work, we will discuss how the 
underlying cosmological model affects the order statistics and con- 
front the theoretically derived order statistics with observations, af- 
terwards. 



3 DEPENDENCE OF THE ORDER STATISTICS ON THE 
UNDERLYING COSMOLOGY 

Eventually, the order statistics in mass and redshift is determined 
by the number of galaxy clusters in a given cosmic volume. The 
quantities that impact on this number can be categorised into two 
classes. The first one contains all effects that modify structure for- 
mation itself, like the choice of the mass function or the amplitude 
of the mass fluctuations, erg, for instance. These effects manifest 
themselves most strongly in the exponentially suppressed tail of the 
mass function, hence at high masses. The second class contains all 
the effects that modify the geometric evolution of the Universe. By 
changing the evolution of the cosmic volume, the number of clus- 
ters in a given redshift range can be substantially different, even if 
both cosmologies yield the same the number density of objects of 
a given mass (see e.g. Pace et al. 2010). 

3.1 Impact of cosmological parameters 

In order to quantify the impact of different cosmological parame- 
ters on the order statistics in mass and redshift, we study the effect 
on the 98-percentile, g98, which we use to define possible outliers 
from the underlying distribution. In Fig. 2, we present the relative 
difference in g98 as a function of four different cosmological pa- 
rameters comprising erg, (assuming the flatness constraint), the 
equation of state parameter, wq, and the derivative Wa from the re- 
lation w{a) - wq + Wa{\ - a), where a denotes the scale factor. A 
non- vanishing value of the latter indicates a time- varying equation 
of state. In each panel of Fig. 2, we show the relative differences 
for 5 different orders, of the order statistic in mass with z e [0, oo] 
(blue lines) and z e [1, oo] (green lines), as well as in redshift (red 
lines) assuming mum = 10^^ M©. For all calculations we assumed 
the full sky and the Tinker et al. (2008) mass function. 

It can be seen that order statistics is very sensitive to erg, such 
that the relative differences in 298 would amount to ~ 7 per cent 
for the range allowed by WMAP7 of (erg = 0.811 + 0.023). All 
three order statistics exhibit the same functional behaviour, with 
the mass-based ones being more sensitive than the redshift-based 
one. This can be understood by the fact that the mass-based order 
statistics probe the most massive clusters and hence the exponential 
tail of the mass function which is highly sensitive to erg. 

For modifications of the matter density, Qm, assuming the flat- 
ness constraint = 1 - fl^, the situation is substantially different 
from the previous case (see upper right panel of Fig. 2). Overall, the 
order statistics are less sensitive and they do not exhibit the same 
functional behaviour. The order statistics in mass (blue lines) per- 
forms best for larger value of flm because the most massive clusters 
will reside at rather low redshifts. At high redshifts (green and red 
lines), the increase in Q.^ and hence, the decrease in Qa, yields a 
smaller number of very massive clusters. Despite the increase in the 
matter density, the decrease in volume is dominating for the range 



of Qm shown in the plot and, thus, the relative difference decreases. 
In this sense the volume effects dominate at high redshifts over the 
increase in matter density, whereas at low redshifts the increase in 
matter density dominates. 

The lower left panel of Fig. 2 shows the sensitivity of the order 
statistics to changes in the constant equation of state, wq. Evidently, 
the most massive clusters at low redshifts (blue line) have no sensi- 
tivity to Wo, whereas at high redshifts (green and red line) the sen- 
sitivity is better. The volume effects are, compared to modifications 
in ^^m, less important and the observed increase in the relative dif- 
ference in g98 with decreasing wq is dominated by modifications 
of the exponential tail of the mass function (for a more thorough 
discussion, see e.g. Pace et al. 2010). 

When assuming a time-dependent equation of state, modelled 
by w{a) - wq + Wa{\ - a), as presented in the lower right panel 
of Fig. 2, the observed functional behaviour can be explained by 
identical arguments as before. The results exhibit again the high 
sensitivity of the high redshift order statistics on modifications of 
Wa. It should be noted that we fixed wq - -1.0 for all cases. 

It can be summarised that for modifications that strongly af- 
fect the structure formation, like erg for instance, the order statistics 
in mass for z e [0, oo[ is comparable in its sensitivity to the redshift 
based order statistics. Modifications that strongly alter the geomet- 
ric evolution of the Universe affect more strongly the order statis- 
tics in redshift. However, one should keep in mind that in the case 
of the order statistics in mass, the relative differences are on the 
same level as the inaccuracies in cluster mass estimates. This prob- 
lem does not occur for redshifts, which can be measured to a very 
high accuracy. Of course, in this case the observational challenge is 
transferred to compiling a sample with a precise mass limit. Apart 
from the cosmological parameters also the choice of the mass func- 
tion is expected to have a strong effect on the order statistics as will 
be discussed in the following subsection. 



3.2 Impact of the choice of the mass function 

When performing a falsification experiment of ACDM using the n 
most massive or n highest redshift clusters, then one has to spec- 
ify the reference model against which the observations have to be 
compared with. Apart from the cosmological parameters that are 
usually fixed to the obvious choice of the WMAP7 values, a halo 
mass function has to be chosen as well. As mentioned earlier, this 
is particularly important for galaxy clusters since the exponentially 
suppressed tail of the mass function is naturally very sensitive to 
modifications. 

In order to quantify the impact of different mass functions 
on the order statistics in mass and redshift, we computed the 
cdfs, F^n-9).---.F{n). for the Press & Schechter (1974) (PS), the 
Tinker et al. (2008) and the Sheth & Tormen (1999) (ST) mass 
functions for /sky = 1 and present them from top to bottom in Fig. 3 . 
Comparing the panels to each other reveals the tremendous sensi- 
tivity of the distributions to the choice of the mass function. Taking 
the Tinker mass function as a reference, the median, 250, changes 
for both types of order statistics by -20 per cent for the PS case 
and by +15 percent for the ST case. These differences can be ex- 
plained by the fact that the ST mass function leads to a substantial 
increase in the number of haloes, particularly at the high mass end, 
whereas the PS mass function results in much fewer haloes in the 
mass and redshift range of interest. For the remainder of this pa- 
per we will use the Tinker mass function as reference because the 
halo masses are defined as spherical overdensities with respect to 
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Figure 2. Relative differences in the 98 percentile, ^98, with respect to the ACDM case for the order statistics in mass (blue lines), in mass with a lower 
redshift limit of z = 1 (green lines) and in redshift (red lines) as a function of different cosmological parameters. The upper left panel shows the variation 
with erg, the upper right panel the one with the lower left one the one with the constant equation of state parameter wq and the lower right one shows the 
variation with the derivative of a linearised model for a time dependent equation of state Wa . The different line-styles denote different orders as indicated in the 
individual panels. For all calculations, the full sky and the Tinker mass function were assumed. 



the mean background density, a definition that is closer to theory 
and actual observations than friend-of-friend masses. 

However, considering that due to statistical limitations, cur- 
rent fits for the mass function are still not very accurate for the 
highest masses (> 3 x 10^^ M©) and that systematic uncertainties 
allow even smaller masses an accuracy of a few per cent at most 
(Bhattacharya et al. 2011), one has to be very cautious with falsi- 
fication experiments that are based on extreme objects. The uncer- 
tainty in the mass function alone will allow a rather wide range of 
distributions. 



4 SUITABLE SAMPLES OF GALAXY CLUSTERS FOR 
AN ORDER STATISTICAL ANALYSIS 

Having introduced the order statistics of the most massive or the 
highest redshift clusters, we intend now to compare observed clus- 
ters with the theoretical distributions. To do so, it is necessary to 
select suitable samples of galaxy clusters, which we will discuss in 
the following. 



4.1 General considerations 

The selection of a suitable sample of galaxy clusters for an or- 
der statistical analysis is by no means a trivial task. The neces- 
sary ordering of the quantities mass and redshift by magnitude 
requires that they have been derived in an identical way across 
the sample. Otherwise, systematics and biases, like the differences 
between lensing and X-ray mass estimates for instance (see e.g. 
Mahdavi et al. 2008; Zhang et al. 2010; Planck Collaboration et al. 
2012; Meneghetti et al. 2010; Rasia et al. 2012), will render the or- 
dering meaningless. Despite an increasing amount of data from 
diff'erent surveys, a lack of large homogeneous samples persists. 
Thus, we decided to base our comparison on clusters that stem from 
catalogues like the SPT massive cluster sample (Williamson et al. 
2011) and the MCXC cluster catalogue (Piff^aretti et al. 2011), 
which will be discussed in further detail below. 



4.2 The SPT massive cluster sample 

The SPT survey (Carlstrom et al. 2011) is ideally suited for the in- 
tended purpose of an order statistical analysis. Being based on the 
Sunyaev Zeldovich (SZ) eff^ect (Sunyaev & Zeldovich 1972, 1980) 
the SPT survey is able to detect massive galaxy clusters up to 
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Figure 3. Impact of the different mass function on the cdfs of the first ten orders F{n), . . . ,F(^n-9) in mass (left column) and redshift (right column). All 
distributions were computed for the full sky and three different mass functions comprising the Press-Schechter (PS), the Tinker and the Sheth& Tormen (ST) 
ones, ordered from top to bottom. For the distributions in redshift, a limiting survey mass of mum = 10^^ M© has been assumed. 



high redshifts. The fact that the limiting mass of SZ surveys varies 
weakly with redshift (Carlstrom et al. 2002) allows in principle to 
construct mass limited cluster catalogues. However, it should be 
emphasised that the assumption of an independent of redshift 
depends critically on the sensitivity and the beam width of an actual 
survey. 

For this work, we take the catalogue of Williamson et al. 
(2011) which comprises the 26 most significant detections in the 
full survey area of Af^^ = 2500 deg^. Ensuring a constant mass 
limit of M2oom =^10^^ M©, clusters were selected on the basis of a 
signal-to-noise (S/N) threshold in the filtered SPT maps. For all 26 
catalogue members, either photometric or spectroscopic redshifts 
were determined as well. The cluster masses given in the catalogue 



are defined with respect to the mean cosmic background density 
and need no further conversion to match the mass definition of the 
reference Tinker et al. (2008) mass function. To each cluster of the 
sample we assign the error bars that we obtained by adding the re- 
ported statistical and systematic errors in quadrature. 

4.3 The MCXC duster catalogue 

The MCXC catalogue (Piffaretti et al. 2011) is based on the pub- 
licly available compilation of clusters' detections from ROSAT 
All-Sky Survey (NORAS, REFLEX, BCS, SGP, NEP, MACS, and 
CIZA) and other serendipitous surveys (160SD, 400SD, SHARC, 
WARPS, and EMSS), and provides the physical properties of 



© 2012 RAS, MNRAS 000, 1-17 



Order statistics of galaxy clusters 1 



Table 1. Compilation of the ten most massive galaxy clusters from the SPT massive cluster sample (Williamson et al. 2011) and the MCXC catalogue 
(Piffaretti et al. 2011), respectively. The masses M2oom and ^^200111 respect to the mean background density before and after the correction for the 



Eddington bias based on the estimated mass uncertainty cr\aM- 
is based on. 



The last column lists the references for the values of the observed mass, on which the analysis 



Rank 


Cluster 


z 


A/f'^r\r\ in unit*; r\f A/f^ 




yi^Edd • of yi/^ 

200m I^AIALS Ui IVlQ 




SPT catalogue (A^ = 2500 de^ 












pt 


SPT-CL J0658-5556 


0.296 


(3.12 ± 1.15) X 10^^ 


0.39 


1.99:0^x1015 


Williamson et al. (2011) 


2nd 


SPT-CL J2248-4431 


0.348 


(2.90+ 1.03) X 10^5 


0.37 


l-9l!n-59 X 10^5 


" 


3rd 


SPT-CL JO 102-49 15 


0.870 


(2.16 ± 0.32) X 10^5 


0.15 


1.98+0|2 X 1015 


Menanteau et al. (2012) 


4th 


SPT-CL J0549-6204 


0.320 


(1.99 + 0.67) X 10^5 


0.35 


1.48:of X 1015 


Williamson et al. (2011) 


5* 


SPT-CL J0638-5358 


0.222 


(1.91 ±0.62) X 10^5 


0.34 


1.50+^-61 X 1015 


" 


6* 


SPT-CL J0232-4421 


0.284 


(1.88 ± 0.59) X 10^5 


0.32 


1.48+^-5^ X IO15 


" 


yth 


SPT-CL J0645-5413 


0.167 


(1.81 +0.60) X 10^5 


0.34 


1.43+|58 X 1015 


" 


gth 


SPT-CL J0245-5302 


0.098 


(1.70 + 0.46) X 10^^ 


0.25 


1.48+11 X 1015 


" 


9th 


SPT-CL J220 1-5956 


0.300 


(1.70 ± 0.42) X 10 


0.28 


1 A 0+0 48 V ^ 1 r\l 5 

1.48 X 10^^ 




10* 


SPT-CL J2344-4243 


0.450 


(1.65 ± 0.38) X 10^5 


0.31 


1.28+0-g X 1015 




MCXC catalogue (A^ = 27490 deg^) 










pt 


J0417.5-1154 


0.4430 


(3.86 + 0.62) X 10^5 


0.15 


3.55:0-57 X 1015 


Piffaretti et al. (2011) 


2nd 


J22 11.7-0349 


0.3970 


(3.21 ±0.51) X 10^5 


0.15 


2.98+0-48 X 1015 




3rd 


J2243.3-0935 


0.4470 


(3.04 ± 0.49) X 10^^ 


0.15 


2.82+035 X 1015 




4th 


J0308.9-F2645 


0.3560 


(2.95 ± 0.47) X 10^^ 


0.15 


2.75+0^45 X 1015 
2.25+o| X 1015 




5* 


J1504. 1-0248 


0.2153 


(2.36 ± 0.35) X 10^5 


0.14 




6* 


J1347.5-1144 


0.4516 


(2.30 ± 0.41) X 10^5 


0.16 


2.14+0^37 X 1015 
2.14+o|5 X 1015 
2.06+o:| X 1015 
2.01+0:33 X 1015 
2.03+o:| X 1015 




yth 


J1731.6-F2251 


0.3890 


(2.28 ± 0.37) X 10^5 


0.15 




gth 


J0717.5-F3745 


0.5460 


(2.21 ± 0.35) X 10^5 


0.15 




9th 


J2248.7-4431 


0.3475 


(2.13 ± 0.35) X 10^5 


0.15 




10* 


J1615.7-0608 


0.2030 


(2.14 ± 0.30) X 10^5 


0.13 





1743 galaxy clusters systematically homogenised to an overden- 
sity of 500 (with respect to the cosmic critical density). This meta- 
catalogue is not complete in any sense, but it is constituted by X-ray 
flux-limited samples that ensure that the X-ray brightest objects in 
the nearby (z < 0.3) Universe, and therefore the most massive X- 
ray detected clusters, are all included. 

We have then simply ranked the objects accordingly to their 
estimated M2oom, that is obtained from the tabulated M500C as 



■Ms, 



200Q, (R 



500 



(13) 



where O.^ = + zfjEl, E, = (^^^(1 + zf + ^aY^^, and the 

ratio between the radii at different overdensities has been obtained 
by assuming an NFW profile (Navarro et al. 1996) with C200 - 4. 



4.4 Preparations of the ordered samples 

We order the SPT and MCXC catalogues by magnitude of the ob- 
served mass and present the ten most massive systems in Table 1 . 
For statistical comparisons the observed masses have to be cor- 
rected for the Eddington bias (Eddington 1913) in mass. As a re- 
sult of the exponentially suppressed tail of the mass function and 
the substantial uncertainties in the mass determination of galaxy 
clusters, it is more likely that lower mass systems scatter up while 
higher mass systems scatter down, resulting in a systematic shift. 
Thus, before an observed mass can be compared to a theoretical 
distribution, this shift has to be corrected for. To do so, we follow 
Mortonson et al. (2011) and shift the observed masses. Mobs, to the 



corrected masses, Mcon, according to 



InMcc 



1 2 

: In Mobs + 2^^1nM' 



where e is the local slope of the mass function (dw/dlnM oc M^) 
and cTinM is the uncertainty in the mass measurement. We corrected 
the observed masses in both, the SPT and the MCXC catalogues, 
using the values of crinM hsted in the fifth column of Table 1 which 
we deduced from the reported uncertainties in the nominal masses. 
The larger the observational errors are, the larger is the correction 
towards lower masses. 

As an exemplary exception from the SPT catalogue, we 
used for the mass of SPT-CL JO 102-49 15 the value reported by 
Menanteau et al. (2012), which is based on a combined SZ-hX- 
rays-hoptical-hinfrared analysis. The multi-wavelength study shifts 
Mobs = (1.89 + 0.45) X 1015 Mo (Williamson et al. 2011) to a larger 
value of Mobs = (2. 16 + 0.32) x 10i5 M©, changing the rank from the 
fifth to the third most massive. This shows that with the expected 
increase in the quality of cluster mass estimates, the ordering of the 
most massive cluster will undergo significant changes. We expect 
that the reshuffling will affect more strongly the most massive clus- 
ters due to the fact that the large error bars will cause lower ranked 
clusters to scatter up. We will discuss the impact of the reshuffling 
in more detail in Sect. 5.1. 

In addition, we sorted the SPT catalogue by redshift and list 
the ten highest redshift clusters above mum ^ 10i5 M© in Table 2. 



(14) 



5 COMPARISON OF THE INDIVIDUAL ORDER 
STATISTICS WITH OBSERVATIONS 

In this section we will compare the individual ranked systems listed 
in Table 1 for the mass and in Table 2 for the redshift with the indi- 
vidual distributions for each rank, as e.g. shown in Fig. 3. 
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Figure 4. Functional box plots for the first nine orders in the observable mass as indicated in the individual panels. Here, the red line denotes the median 
(250), the blue-bordered, grey region the interquartile range (IQR) and the black lines denote the 2 and the 98-percentile (^2, ^98). The green error bars 
show the corresponding observed masses, M200m, from the SPT (green circles) and the MCXC (green triangles) catalogues (see Table 1) for their respective 
survey areas of Af^ = 2500 deg^ and Af^^^ = 27490 deg^. 



Table 2. Compilation of the ten highest redshift clusters from the SPT mas- 
sive cluster sample (Williamson et al. 2011). Here, (s) and (p) denote the 
spectroscopic and photometric redshifts, respectively. 



Rank 


Cluster 


z 


^200m in units of M© 




SPT-CL J2106-5844 


1.132 (s) 


(1.27 ± 0.21) X 


10^5 


2nd 


SPT-CL J0615-5746 


0.972 (s) 


(1.32 ± 0.40) X 


1015 


3rd 


SPT-CL JO 102-49 15 


0.870 (s) 


(2.16 ± 0.32) X 


1015 


4th 


SPT-CL J2337-5942 


0.775 (s) 


(1.99 ± 0.20) X 


1015 


5* 


SPT-CL J2344-4243 


0.620 (p) 


(1.91 +0.50) X 


1015 


6* 


SPT-CL J0417-4748 


0.620 (p) 


(1.88 + 0.20) X 


1015 


yth 


SPT-CL J0243-4833 


0.530 (p) 


(1.81 ±0.23) X 


1015 


gth 


SPT-CL J0304-4401 


0.520 (p) 


(1.70 ± 0.33) X 


1015 


gth 


SPT-CLJ0438-5419 


0.450 (p) 


(1.70 + 0.38) X 


1015 


10* 


SPT-CLJ0254-5856 


0.438 (s) 


(1.65 ± 0.25) X 


1015 



5.1 Order statistics in cluster mass 

In order to demonstrate the impact of the survey area on the dis- 
tributions of the order statistics in mass, we show in Fig. 4 the 
dependence of different quantiles (22, 225, 250, 275 and 298) 
on the survey area for the nine most massive clusters. In addi- 
tion, the green error bars show the clusters from the SPT and 
MCXC catalogues listed in Table 1 for the respective survey areas 
of AfT = 2500 deg2 and Af^^^ = 27490 deg^. 

From the individual panels in Fig. 4 it can be inferred that, 
as expected, a larger survey area yields a larger expected mass 
for the individual rank. Furthermore, with increasing rank towards 
higher orders, the interquantile range, like (Q2-Q98), narrows. A 
behaviour that can also be seen in Fig. 1 as steepening of the cdf 
with increasing rank. Therefore, the largest mass (first order) is ex- 
pected to be realised in a much wider mass range than the higher 
orders. 
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Figure 5. Box-and- whisker diagram of the ten most massive clusters from the SPT survey (left column) and the MCXC catalogue (right column) for three 
different choices of the mass function as denoted in the title of each panel. For each order, the red lines denote the median (^50), the blue-bordered, grey 
boxes give the IQR and the black whiskers mark the range between the 2 and 98-percentile (Q2, ^98) of the theoretical distribution. The green, filled circles 
denote the nominal observed cluster masses, M200m, the orange, empty triangles the ones that are corrected for the Eddington bias in mass and the violet, 
empty circles are the results of the Monte Carlo reshuffling of the ranks. All error bars denote the Icr range. 



We will now compare the observations in more detail with 
the theoretical expectations in the form of box-and-whisker dia- 
grams as shown in Fig. 5. Here, the blue-bordered, grey filled box 
denotes the interquartile range (IQR) which is bounded by the 25 
and 75 -percentiles (225, 275) and the median (250) is depicted 
as a red line. The black whiskers denote the 2 and 98-percentiles 
(22,298) and we follow the convention that observations that fall 
outside are considered as outliers. As before, the nominal observed 
cluster masses are denoted as green error bars where for the left 
column the SPT catalogue and for the right column the MCXC cat- 
alogue was used. In addition we plot the Eddington bias corrected 



masses, M2qq^, from the sixth column of Table 1 as orange triangles 
with dashed error bars. We performed the analysis for three differ- 
ent mass functions, comprising from the top to the bottom panel, 
the PS, the Tinker and the ST mass functions. In addition to the 
Eddington bias in mass, we expect a shift to larger masses caused 
by the reshuffling of orders due to the uncertainties in mass. In or- 
der to quantify this effect, we Monte Carlo (MC) simulated 10000 
realisations of the 26 SPT and 123 MCXC (with M > lO^^ M©) 
cluster masses after their correction for the Eddington bias and or- 
der them by mass. The masses were randomly drawn from the indi- 
vidual error interval, assuming Gaussian distributions. We present 
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turn, by assuming the ACDM reference cosmology, the compari- 
son of the observed masses with the individual order distributions 
allows to check the completeness of the observed sample. 

In the upper panel of Fig. 6, we present the dependence of 
different percentiles (22,225,250,275 and 298) on the order for 
a survey area of As = 20000deg^. Choosing the 298 percentile 
as exclusion criterion, one would need roughly to find ten clusters 
with m > 2.5 X 10^^ M©, three clusters with m > 3.2 x 10^^ M© or 
one cluster with m > 5 x 10^^ M© in order to report a significant 
deviation from the ACDM expectations. Of course, the observed 
masses might have to be corrected for the Eddington bias in mass 
and a possible reshuffling as previously demonstrated. In general, 
exclusion criteria based on order statistics extend previous works 
(Mortonson et al. 2011; Waizmann et al. 2012a) from statements 
about single objects to statements about object samples which con- 
siderably improves the reliability of the entire study. 



Figure 6. Dependence of different percentiles on the order for the order 
statistics in mass (upper panel) and redshift (lower panel). For both cases 
a survey area of As = 20000deg^ is assumed. Further, a limiting mass 
of miim = 10^^ Mq has been adopted for the order statistics in redshift. 
All percentiles are denoted by the same line styles as used in the previous 
figures and are given in the key. 

the results as violet, empty circles with dash-dotted Icr error bars 
in Fig. 5. It can be seen that the highest ranks are more strongly 
affected by the reshuffling than the lower ones and that they are 
on average shifted to larger values. Of course, the amount of this 
effect will depend on the size of the error bars. Further, the reshuf- 
fling yields mass values that fall between the nominal (green error 
bars) and the Eddington bias corrected ones (orange error bars). 

For the SPT catalogue, it can be seen from the top left panel of 
Fig. 5 that the outdated PS mass function seems to be disfavoured 
by the reshuffled and the nominal masses of the ten largest objects. 
However, the error bars are large and do not allow an exclusion of 
the PS mass function. For the Tinker and the ST mass function, 
the boxes indicating the theoretical distributions move to larger 
mass values and therefore they match the observed masses better 
than the PS mass function. In particular, the third ranked (second 
ranked after Eddington bias correction) system SPT-CL JO 102 with 
its smaller errors and, hence, giving the tightest constraints, is con- 
sistent with ACDM for both mass functions. All other ranks are 
consistent as well due to their large error bars. The reshuffled sam- 
ple matches perfectly the Tinker mass function consolidating the 
conclusion that the most massive clusters of the SPT sample are 
in agreement with the statistical expectations. The conclusions for 
the MCXC catalogue are identical, however the jump between the 
fourth and the fifth largest order yield to an inconsistency of the 
observed higher orders with the expectations based on the ST mass 
function. This jump is clearly caused by the incompleteness of the 
MCXC catalogue and, thus, the inclusion of the missing clusters 
would most certainly move the observed sample to higher masses 
in the direction of the results we obtained from the SPT sample. In 
this sense we do not see any indication of a substantial difference 
between the small and wide field survey. 

The analysis of the SPT sample illustrates the potential of util- 
ising the n most massive galaxy clusters to test underlying assump- 
tions, like e.g. the mass function. For instance, a multi-wavelength 
study of the 26 SPT clusters would reduce the error bars to the level 
of SPT-CL JO 102 (the nominal third ranked cluster in the left col- 
umn of Fig. 5 ), which would significantly tighten the constraints 
on the underlying assumptions like e.g. the halo mass function. In 



5.2 Order statistics in duster redshift 

We performed an identical analysis for the individual order statis- 
tics for the SPT massive cluster catalogue ranked by redshift listed 
in Table 2. For the theoretical distributions we assume a limiting 
mass of miim = 10^^ M© and a survey area of A^^^ - 2500 deg^. As 
before, we present in Fig. 7 the dependence of the order statistical 
distributions on the survey area for the first nine orders. Again, an 
increase in the survey area yields a shift of the theoretical distribu- 
tions to higher redshifts and, as shown in the right panel of Fig. 1 , 
the cdfs steepen for the higher ranks, resulting in a shrinking in- 
terquantile range. 

In Fig. 8, we present the box-and- whisker diagram in redshift, 
again for the PS, the Tinker and the ST mass functions (from top to 
bottom). The definition of boxes and whiskers remains unchanged 
with respect to Fig. 5. Again, the data from Table 2 is denoted by 
green error bars, which are negligibly small in the case of spectro- 
scopic redshifts. Thus, we abstained from the MC simulation of the 
reshuffling in the case of redshift. While for the order statistics in 
mass the results only depended on the choice of the survey area, 
the situation is different for the order statistics in redshift. Here, 
a constant survey limiting mass is assumed, which will be subject 
to uncertainties for a real survey and, furthermore, will also exhibit 
some redshift dependence. Thus, the theoretical distributions are in- 
trinsically less accurate than the ones with respect to cluster mass. 
Indeed, the comparison with the data in Fig. 8 exhibits a different 
behaviour with respect to the one in Fig. 5. Here, first four orders 
seem to be fit better by the Tinker mass function while the higher 
orders seem to favour the PS mass function. Taking the Tinker mass 
function as reference it seems that a few systems with M > 10^^ M© 
are missing at redshifts z > 0.7. The difference with respect to the 
findings for the order statistics in mass for the same sample could, 
along the lines of Sect. 3, be interpreted as a signature of a devi- 
ation from the reference ACDM model. However, considering the 
previously mentioned simplifying assumptions in the modelling of 
the theoretical distributions, we do not infer any cosmological con- 
clusions and leave a better, more realistic, modelling of miim(z) of 
the SPT survey to a future work. 

In the lower panel of Fig. 6, we present the dependence of 
different percentiles (22,225,250,275 and 298) on the order for 
a survey area of As = 20 000 deg^ and a constant limiting mass of 
^lim = 10^^ Mq. Taking the 298 percentile as exclusion criterion, 
one would need to find ten clusters with z > 1, three clusters with 
z> 1.2 or one cluster with z> 1.55 in order to report a significant 
deviation from the ACDM expectations. Currently, SPT-CL J2106 
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Figure 7. Functional box plot for the first nine orders in the observable redshift as indicated in the individual panels. Here, the red line denotes the median 
((250), the blue-bordered, grey region the interquartile range (IQR) and the black lines denote the 2 and the 98-percentile (^2, ^98). The green circles (with 
error bars in the case of photometric redshifts) denote the redshifts from the SPT catalogue as listed in Table 2. 



is the only known cluster of such a high mass having a redshift z > 
1. With an assigned survey area of A, = 2800 deg^ (ACT+SPT), 
it might from a statistical point of view still be possible to find ten 
objects that massive at z > 1 in the larger survey area. The method 
presented in this work allows to construct similar exclusion criteria 
for any kind of survey design. 



6 COMPARISON OF THE JOINT ORDER STATISTICS 
WITH OBSERVATIONS 

Having studied the individual order statistics in mass and redshift 
in the previous section, we turn now to the study of the joint distri- 
butions of the order statistics as introduced in Sect. 2. 1 . 

The simplest case of a joint order distribution is two- 
dimensional. In this case the pdf and cdf are given by equation 6 
and equation 7, respectively. Starting with the joint pdf, we present 
in Fig. 9 the joint distributions in mass (left panel) and redshift 



(right panel) for several order combinations as denoted in the in- 
dividual panels. All calculations assume the full sky and the Tinker 
mass function. In the case of the joint distributions in redshift, we 
assume a constant limiting survey mass of mum = 10^^ M©. Due to 
the condition that x <y, all distributions are limited to a triangular 
domain. 

An inspection of the diff'erent pdfs in Fig. 9 reveals that, for 
the combination of the first and the second largest order (upper left- 
most panel), the most likely combination of the observables is very 
close to the diagonal. This means that it is more likely to find the 
two largest values close to each other, at absolute values that are 
smaller than the extreme value statistics would imply for the max- 
imum alone. Then, when moving to combinations of the first with 
higher orders (first row), it can be seen that the peaks of the pdfs 
move away from the diagonal and that they extend to larger values 
for the larger observable. This indicates that it is more likely to find 
the two systems with a larger separation in the observable when the 
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Figure 8. Box-and- whisker diagram of the ten highest redshift clusters from 
the SPT catalogue as listed in Table 2 for three different choices of the mass 
function as denoted in the title of each panel. For each order, the red lines 
denote the median (^50), the blue-bordered, grey boxes give the IQR and 
the black whiskers mark the range between the 2 and 98-percentile 
298). The green circles (with error bars in the case of photometric redshifts) 
represent the observed cluster redshifts. 



difference between the considered orders is larger. Accordingly, for 
higher order combinations (lower rows), the peaks of the joint pdf 
move to smaller values of the observables. It should also be noted 
that the peaks steepen for higher order combinations, confining the 
pdfs to smaller and smaller regions in the observable plane. As an 
example, the first and second largest observations (upper leftmost 
panel) can be realised in much larger area than the sixth and eighth 
largest one (lower rightmost panel). 

Apart from the joint pdfs, it is also instructive to study the joint 



cdfs as presented in Fig. 10 for the observed mass (left panel) and 
redshift (right panel). In order to add observational data from the 
SPT catalogue, we assume a survey area of = 2500 deg^ and 
a miim = 10^^ Mq for the joint distribution in redshift. Additionally, 
we added the two largest nominal observed (red error bars) and 
the Eddington bias corrected masses (grey error bars) from Table 1 
to the left panel and the two highest redshifts of the SPT massive 
cluster sample from Table 2 to the left panel. In the case of the 
mass, we find F(^n-\){n) ~ 0.92 for the nominal and F(^n-i){n) ~ 0.1 
for the Eddington bias corrected masses. Hence, using the central 
values, in ~ (8 - 90) percent of the cases a mass larger than the one 
of SPT-CL J0658 and a mass larger than the one of SPT-CL J2248 
are observed. Thus, also the joint cdf confirms that the two largest 
masses do not exhibit any tension with the concordance cosmology. 
The same conclusion applies in the case of the joint distribution in 
redshift. 

By means of equation 8 these results can be extended to the n- 
dimensional case, allowing the formulation of a likelihood function 
of the ordered sample of the n most massive or highest redshift 
clusters. 



7 SUMMARY 

In this work, we studied the application of order statistics to the 
mass and redshifts of galaxy clusters and compared the theoreti- 
cally derived distributions with observed samples of galaxy clus- 
ters. Our work extends previous studies that hitherto considered 
only the extreme value distributions in mass or redshift. 

On the theoretical side, our results can be summarised as fol- 
lows. 

(i) We introduce all relations necessary to calculate pdfs and 
cdfs of the individual and joint order statistics in mass and red- 
shift. In particular, we find a steepening of the cdfs for higher or- 
der distributions with respect to the extreme value distribution of 
both mass and redshift. This steepening corresponds to a higher 
constraining power from distributions of the ^-largest observations. 
The presented method extends previous works to include exclusion 
criteria based on the n most massive or n highest redshift clusters 
for a given survey set-up. 

(ii) Conceptually, we avoid the bias due to an a posteriori choice 
of the redshift interval in the case of the order statistics in mass by 
selecting the interval < z < 00. Hence, we study the statistics 
of the hierarchy of the most massive haloes in the Universe, which 
mostly stem from redshifts z < 0.5. On the contrary, when choos- 
ing the order statistics in redshift, focus is laid on haloes that stem 
from the highest possible redshifts. However, the calculations will 
require a model of the survey characteristics in the form of a limit- 
ing survey mass as a function of redshift. 

(iii) By putting the emphasis on either the most massive or on 
the highest redshift clusters above a given mass limit, the order 
statistics is e.g. particularly sensitive to the choice of the mass func- 
tion. While the order statistics in mass is very sensitive to erg and 
Qm due to the domination of low redshift objects, the order statistics 
in redshift proves to be very sensitive to wq and Wa. For a fixed cos- 
mology, both order statistics are efficient probes of the functional 
shape of the mass function at the high mass end. 

(iv) In addition to the individual order statistics, we study as ex- 
ample case also the joint two order statistics. We find that for the 
combination of the largest and the second largest observation, it is 
most likely to find them to be realised with very similar values with 
a relatively broadly peaked distribution. 
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the joint distribution in redshift. The color bar is set to range from to the maximum of the joint pdf for each rank combination. 
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Figure 10. Joint cdf F(^)(^)(x, j) (see equation 7) for the observable mass (left panel) and redshift (right panel) for the combination of the largest and second 
largest observation, assuming a survey area of - 2500 deg^ and a constant limiting survey mass of mum - 10^^ for the joint distribution in redshift. 
The red and grey error bars denote the nominal and the Eddington bias corrected values, M200m and >W200m' catalogue as listed in Table 1 for the 

mass and in Table 2 for the redshift. 



(v) In order to allow a quick estimation of the distributions of 
the order statistics, we provide in the Appendix B fitting formulae 
for F{x) as a function of the survey area for the cases of mass and 
redshift. The fitting formulae for the distribution in mass allow for a 
percentile estimation in the range from Ql to g98 with an accuracy 
better than one per cent for As > 200 deg^ and for the ten largest 
masses. In the case of the order statistics in redshift, the quality of 
the fits depends on the chosen mum. However, for survey areas of 
As > 2000 deg^ accuracies better than two per cent can be achieved 
for large values of mum = 10^^ M© and a lowering of mum further 
improves the accuracy. 

After introducing the theoretical framework, we compared the the- 
oretical distributions with actually observed samples of galaxy 



clusters that we ranked by the magnitude of the observables mass 
and redshift. We decided to compile two catalogues, the main one is 
based on the SPT massive cluster sample (Williamson et al. 2011) 
and additionally we analysed the meta-catalogue of X-ray detected 
clusters of galaxies MCXC (Piff'aretti et al. 2011) based on pub- 
licly available flux-limited all- sky survey and serendipitous cluster 
catalogues. This meta-catalogue can be considered as complete for 
z < 0.3 and, hence, by no means as complete as the SPT one. The 
results of the comparison can be summarised as follows. 

(i) In the case of the order statistics in mass, we compared the 
theoretical expectations for the ten largest masses for the PS, the 
Tinker and the ST mass functions. Assuming WMAP7 parameters, 
we find that the nominal and the Eddington bias corrected values 
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for the observed masses favour the Tinker and the ST mass func- 
tions. When considering the possible bias due to a reshuffling of the 
ranks caused by the large error bars (statistical + systematic errors), 
we find that the SPT sample matches the Tinker mass function very 
well. The constraints are expected to tighten considerably once the 
error bars of all objects are scaled down by combining several clus- 
ter observables in multi- wavelength studies. 

(ii) In contrast to the ranking in mass, the order statistics of the 
SPT clusters in redshift is less well fit by the theoretical distribu- 
tions based on the Tinker mass function. It appears that a few sys- 
tems with M > 10^^ Mq are missing at redshifts z > 0.7. One expla- 
nation could be found in a non-standard cosmological evolution to 
which the order statistics in redshift is more sensitive. However, it 
is more likely that a more precise modelling (including the redshift 
dependence) of the true limiting survey mass of SPT will account 
for the observed deviations. 

(iii) Instead of utilising order statistics to perform exclusion ex- 
periments, it can also be used for consistency checks of the com- 
pleteness of the observed sample and of the modelling of the survey 
selection function as indicated by the analysis of the MCXC (mass) 
and the SPT (redshift) samples. 



8 CONCLUSIONS 

We introduced a powerful theoretical framework which allows to 
calculate the expected individual and joint distribution functions of 
the ^-largest masses or the ^-highest redshifts of galaxy clusters 
in a given survey area. This approach is more powerful than the 
extreme value statistics that focusses on the statistics of the single 
largest observation alone. 

As a proof of concept, we compared the theoretical distribu- 
tions with observed samples of galaxy clusters. However, data of 
sufficient quantity, uniformity and completeness is still sparse such 
that constraints are not particularly tight. This situation will most 
certainly improve in the near and intermediate future. Since the 
emphasis of this work lies on the introduction of the theoretical 
framework of order statistics and its application to galaxy clusters, 
we contended ourselves with a study of cluster masses and red- 
shifts. Unfortunately, the mass of a galaxy cluster is not a direct 
observable and subject to large scatter and observational biases. In 
a follow-up work, we intend to extend the formalism to direct ob- 
servables, like for instance X-ray luminosities, and to include the 
scatter in the scaling relations into the theoretical distributions. 
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APPENDIX A: ORDER STATISTICS 

In this appendix, we outline the derivation of the most impor- 
tant relations of the order statistics and some subtleties consid- 
ering their implementation. For more details we refer to the ex- 
cellent textbooks on the topic by Arnold et al. (1992) and by 
David & Nagaraja (2003) which we closely follow for the remain- 
der of this appendix. 
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Figure Al. Schematic for the derivation of f{i){,x) and /(r)(5)(x,j). 



Al Individual distributions 

Let Xi, X2, . . . , X„ be a random sample of a continuous population 
with the cumulative distribution function, F{x). Further, let < 
X(2) < • • • < X{n) be the order statistic, the random variates ordered 
by magnitude, where is the smallest (minimum) and de- 
notes the largest (maximum) variate. The event x < < x + (5x is 
the same as the one depicted in panel {a) of Fig. Al and, thus, we 
have ^ X for / - 1 of the X^, exactly one in x < ^ x + 6x 
and the remaining n- iof the Xi'mXi> x + 6x. Now, the number 
of ways how n observations can be arranged in the three regimes is 
given by 



A{n, i) ■■ 



(/- l)!l!(^-0! ' 
where each of them has a probability of 

[F(x)r' [F(x + Sx) - F(x)] [1 - F(x)r-' . 



(Al) 



(A2) 



Therefore, under the assumption that 6x is small, we find for the 
probability 

Pr{x < ^x + Sx]= A(n, i) [F(x)r^ [1 - F(x)T-' f(x)Sx, (A3) 

neglecting terms of 0(6 x)^. Dividing by 6x and performing 6x^0 
yields the pdf as given in equation 1 



Mx) = lim^ 



Pr{x < X(i) x + 6x] 



= A(nJ)[F(x)r' [1 - F(x)r-^ fix). 



(A4) 



The corresponding cdf of the i-th order, as given by equation 2 in 
Sect. 2, can now either be obtained by integrating the above equa- 
tion or by the following argument 

F(,)(x) = Pr{X(,) < x] 

- Pr{at least / of . . . , are at most x] 

n 

- ^ Pr{exactly k of X(^n)^ • • • , ^(/) are at most x) 

k=i 

= Yj{^^[F(x)f[l-F(x)r-', (A5) 

for -00 < X < 00. Hence, the cdf of is equivalent to the tail 
probability (starting from /) of a binomial distribution with n trials 
and a success probability of F(x). By setting i - n or i - 1 one 
obtains the cdfs for the smallest and the largest order statistics as 
given by equation 3 and equation 4. 



A2 Joint distributions 

The joint pdf of the two order statistics {I ^ r < s ^ n) 

for X < _y can be derived by similar arguments as for the single 
order statistics. The derivation scheme is now extended according 
to panel {b) of Fig. Al. Analogously to equation A4 we obtain 

I Pr{x < Xi^r) <^ X + 6x, y < < j + ] 



f{r){s){x, y) 



lim \ 

6x^0 



6x6y 
A{n, r, s) 

X [F{x)Y-' [F(y) - F{x)Y-'-' 

xnx)ny). 



[1 - F(y)T 



(A6) 



where 

'"' ^> " 7 KIT '^—^r TT • 

(r - 1)!(^ - r - l)l{n - s)\ 

The joint cumulative distribution function can be obtained by 
integrating the pdf from above or again by the following direct ar- 
gument 



F(r)(s)(x, y) 



Pr{X^r) < X, < y] 

Pr{dX least rX(/) < x A at least < y] 

n j 

y y Prfexactly ^ x A exactly jXt^D < y] 



j=s i=r 



|jU/!(7-0!(^-7)! 
[F{x)]' [F{y) - F{x)r' [1 - F{y)f-^ . 



(A8) 



This is exactly identical to the tail probability of a bivariate binom- 
inal distribution. 

Following the same line of reasoning as for the joint two order 
statistics, the above relations can be generalised to the joint pdf of 
Xny,. . . ,Xn^ {I Hi < • • • < Uk ^ n) for x\ • Xk, which is 
given by 

f(xi)-(xk)(X\,- 



{fii - 1)1(^2 -ni-\)\---{n- nk)l 
X [F(x,)r-' fix,) IF{X2) - F{x,)r-^'-' 
xf{x2)"-[l-F{xk)r-''^f{xk). (A9) 

The right hand side of this relation can be written in a more compact 
form (David & Nagaraja 2003) as 

-\nj+i-nj-l \ 



j=i 



\F(xj^,)-F(xj)\ 



7+1 



(AlO) 



where we defined no = 0, n^+i = n + 1, xo = -oo and x^ = +00. 
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Figure A2. Relative differences, A = {Q^^ - Q^^'')IQ^^\ between the fitted 
and directly calculated percentiles (different line styles) as a function of the 
survey area for the order statistics in mass. The differences are shown for 
three different ranks, the largest (black lines), the fifth largest (blue lines) 
and the tenth largest (green lines) one. 



A3 Regarding the implementation 

The implementation of the order statistics for the intended applica- 
tion of this work, as discussed in Sect. 2.1 and Sect. 2.2, is rather 
straightforward. However, one important subtlety arises from the 
combinatoric prefactors that contain factorials of n, which due to 
the large number of haloes cannot be calculated directly. However, 
for all prefactors the factorials of n can be avoided by writing them 
as products and by dividing out common terms. As a simple exam- 
ple, we take the pref actor from equation Al. In this case the index 
/ will, depending on the order, be given by a term like i - {n - j) 
with J = for the distribution of the maximum, j - \ for the second 
largest and so on. Thus, we obtain 



A(n, i - n - j) 



{i - \)\\\{n - i)\ {n- j-\)\{n-n + j)\ 



(All) 



which can be calculated for rather large values of n. In a similar 
manner, all combinatoric prefactors can be simplified and imple- 
mented. 



APPENDIX B: A FITTING FUNCTION FOR THE ORDER 
STATISTICS 

In this additional section, fitting functions for the order statistics in 
mass and in redshift are defined. As functional form for the numer- 
ical fits, we will use equation 5 in combination with the relation 
equation 4, which yields 

-1/7(3;) V^/"^^^ 

i+r(y) 



exp 



''-a(y) 

Piy) 



(Bl) 



Here, x is the observable, either mass or redshift, and the GEV 
parameters or, p and 7 as well as the number of haloes^ , n, are func- 
tions of the survey area via the variable y - log (As). Once the cdf, 

^ For the numerical calculation of n, we limit the mass range without loss 
of generality to the interval relevant for galaxy clusters of 10^^ M© < m < 
10^6 Mo 



F(x), is known, all order statistics can be calculated by means of the 
relations discussed in the previous Appendix A. Inverting the cdfs 
of order statistics allows to obtain the percentiles which can then 
be utilised as ACDM exclusion criteria (see e.g. Fig. 6). 



Bl Order statistics in mass 

In order to determine the fitting function for the order statistics in 
mass, we calculate the GEV parameters according to Davis et al. 
(201 1) and Waizmann et al. (201 1) and the number of haloes, n, as 
a function of the survey area and fit them by the following functions 



aiy) = 5.99888 ln(/ 



10.5689) , 



piy) = 0.362939 exp(-l. 11069/^2^2^^), 
yiy) = -0.239274 ln(3;-0-4^^oo9 + 0.747006) , 

^(3;)= 10^+2.94112^ 



(B2) 
(B3) 
(B4) 
(B5) 



where y - Xog^^iA^). The observable x in equation Bl is defined to 
be X = \ogiQ{M2oomh). We present the results in Fig. A2 in the form 
of relative diff'erences between the fitted and directly calculated val- 
ues of five selected percentiles (22, 225, 250, 275 and 298) as a 
function of the survey area. The diff'erent colors denote the largest 
order statistics, F(„)(x) (black lines), the fifth largest order statistics, 
F(„_4)(x) (blue lines) and the tenth largest order statistics, F(„_9)(x). 
The relative errors in the five diff'erent percentiles are for almost 
the complete range of survey areas on the sub-per cent level (only 
298 for F(„)(x) exhibits a slightly larger error for very small survey 
areas). 



B2 Order statistics in redshift 

For fitting the order statistic in redshift, we proceed in a similar 
way as for the mass, setting x = z in equation B 1 . For calcu- 
lating the GEV parameters as a function of the survey area, we 
follow the approach presented in (Metcalf & Waizmann in prepa- 
ration). However, since in contrast to the order statistics in mass, 
the distributions depend on the choice of the limiting survey mass, 
we fitted the distributions for two choices of miim. First, we set 
^lim = 10^^ Mq, identical to the setup we discussed in this paper 
for the SPT massive cluster sample. Secondly, we lower the thresh- 
old to miim = 5x10^"^ Mq. In the first case, we obtain 



a{y) = 1.13729 ln(0.567735j + 0.332933) , 

piy) = exp[-exp(-1.76728j-^-^^^32 + 0.929307)] , 

yiy) = -2.23597 ln(};-2-^^^^^ + 1.01017) , 

^(-y) ^ -^q0.981095);-1.52015 ^ 

and for the second choice we find 

aiy) = 2.1084 ln(0.284062}; + 1.09002) , 

piy) = exp[- exp(-0.905364};-°-3^^22^ + 1.30066)] , 

y{y) = -0.260275 Iniy'^-^^^^^ + 1.0592) , 

niy) = 100-9985 523;-0.451364 ^ 



(B6) 
(B7) 
(B8) 
(B9) 

(BIO) 
(Bll) 
(B12) 
(B13) 



where y = log^oC^s) for both cases. We present the results in 
Fig. Bl again as relative diff'erences. It can be seen that in the 
case of high limiting mass (upper panel), the fit performs poorly 
for survey areas smaller than ~ 1000 deg^ due to the insufficient 
number of haloes that are expected to be found. However, above 
~ 2000 deg^ the percentiles of the first ten orders can be fitted with 
an accuracy better than two per cent. 
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Figure Bl. Relative differences, A = {Q^^ - 2 )/e , between the fitted 
and directly calculated percentiles (different line styles) as a function of the 
survey area for the order statistics in redshift assuming = 10^^ M© 
(upper panel) and miim = 5 x 10^^ M© (lower panel). The differences are 
shown for three different ranks, the largest (black lines), the fifth largest 
(blue lines) and the tenth largest (green lines) one. 



If the limiting mass is lowered, the quality of the fit improves 
drastically as shown in the lower panel of Fig. B 1 for mum = 
5 X 10^"^ Mq. In this case, sub-percent-level accuracy is reached for 
As > 1000 deg^ and an accuracy better than two per cent down to 
lOOdegl 
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